First order phase transition in Ising model on two connected Barabasi- Albert 

networks. 
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We investigate the behavior of the Ising model on two connected Barbasi-Albert scale-free net- 
works. We extend previous analysis and show that a first order temperature-driven phase transition 
occurs in such system. The transition between antiparalelly ordered networks to paralelly ordered 
networks is shown to be discontinuous. We calculate the critical temperature. We confirm the 
calculations with numeric simulations using Monte-Carlo methods. 
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Phase transitions are one of most interesting phenom- 
ena. While the behavior of systems in noncritical regions 
may be also interesting, most crucial changes appear in 
critical regions. It is therefore important to know when 
such transitions occur, and how do they occur. 
It is widely known that the classic Ising model dis- 
plays a second order temperature driven phase transi- 
tion. The model is throughly investigated, but only on 
regular lattices. Along with the emergence of complex 
networks science starting with breakthrough Barabasi- 
Albert's paper [D, came study of Ising model in such 
systems @, HI, ULTE H| ■ Many aspects of the model have 
been studied, from simple antiferromagnetic interactions 
and spin-glasses 0, HI to the directed structure of the 
network [9(. 

In our previous work [Ic| we have investigated the model 
on a pair of connected networks. Recent research indi- 
cates that one of two phase transitions in such a system 
is in fact a first order phase transition, not second order 
like was thought before. 

In this paper, we investigate the phase transition in a 
pair of connected networks, show evidence that it is in 
fact of first order, and back up our analytical calculations 
with numerical simulations. 



II. MODEL 

In our study, we consider two interconnected Barabasi- 
Albert (B-A) networks, where at each node we place an 
Ising spin. The interactions between the spins are ferro- 
magnetic only. 

The B-A model is a model of a growing network [l| . To 
obtain such a network, one starts with m fully connected 
nodes, and adds new nodes to the network. Each new 
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FIG. 1: Two connected B-A networks. A few nodes from each 
network are shown. The intra-network degrees Icaa and kss 
as well as inter-network degrees kAB and ksB for two sample 
nodes are presented. 



node creates m connections to the existing network. The 
probability that a connection will be made to a node i is 
proportional to its degree This results in a scale- free 
network, with a degree distribution P{k) ~ k~ 3 . 
Our two B-A networks are interconnected by Eab links 
(FigO}. Each of these links connects a node in network 
A with a node in network B. The nodes to be connected 
are chosen preferentially, i.e. the probability to pick a 
given node i equals II Ai = ^aa% / ^AAj ■ If we per- 
form linking in this way, the inter-network degree Uabi 
of a node is statistically proportional its to intra-network 
degree k A Ai- 



III. PHASE TRANSITIONS 

The problem of the Ising model on coupled B-A net- 
works has been considered before [10(. In connected B- 
A networks, Ising model is characterized by two phase 
transition in two different critical temperatures T c _ and 
T c+ . Below T c _ there are two possible phases: both net- 
works ordered in with same spin and both networks or- 
dered with opposite spins. At critical temperature T c _ 
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the state with antiparallel spin ordering disappears, and 
between T c _ and T c+ the system orders only parallely. At 
T c+ and above the temperature is too high for network 
to remain ordered and it assumes paramagnetic state. 
As in regular Ising model, the transition at T c+ is second 
order phase transition. However, unlike previous research 
indicated [l(| the transition at T c _ turns out to be of first 
order. 

We have performed analytic calculations, numeric map 
iterations and Monte-Carlo simulations. 



IV. ANALYTIC APPROACH 




We use a mean-field approach to the problem of Ising 
model. In such approach the self-consistent equation for 
the average spin 



FIG. 2: Hyperbolic tangent plot. Dashed line is Eq[5] for 
T < .T c , while the solid line is for T = T c . Thin lines are axes 
and y = x line. 
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where (3 = 1/T, the temperature T is measured in units 
of inverse Boltzmann constant averaging is over 

the canonical ensemble and hi is the external field acting 
on node i. 

If we consider two networks that interact, we can treat 
the influence of the second network as external field hi. 
Since the inter-network links number is proportional to 
intra-network links we can write the full set of equations 
for two networks 
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We introduce a weighted average spin S = V-E^i kiSi, 
that is an order parameter for the Ising model on ran- 
dom network of nonhomogenous degree distribution. Ad- 
ditionally we put k A Bi = PAkAi, k B Ai = VBksi using 
the fact that inter-network degrees are proportional to 
intra-network degrees (see Sect|H]). We obtain following 
equations for the weighted average spins 

Sa = ^ ~]j<2 tanh (P J AAkAiSA + PJBAkABiS B ) , (5) 



Sb = -pr- tanh {(HBBkBiSB + PJABksAiSA) ■ (6) 



In such a system, a possibility of a first order phase 
transition exists. 

Let us consider a pair of random networks of the same 
size, the same link density and k — const. The right 
side of the Equation [5] is hyperbolic tangent, shifted by 
the value H — JIzapaSb along the x axis. When the 
temperature T is low, (3 is high and the tangent has 
three solutions. If temperature increases to critical T c , 
the curve becomes tangential to the y — x line (Fig(2]). 
The value of H decreases, since network B is also less 
ordered at higher temperature so its influence decreases. 
Below and at T c we can tell that Sb = —Sa from the 
symmetry of the system. At T c , the system is unstable 
and minimal fluctuation of either Sa or Sb causes 
system to switch over to parallel state. 

At T c , the tangent is tangential to the y — x line. We 
can write the conditions for [3 C and Sac 



tanh(f3 c Jk A SA + !3 c Jk A PASB) 
Sa 

dt&nh(f3 c JkASA + PcJIcapaSb) 
dS2 
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We have calculated the (3 C and Sac from these equations 
for Sb = —Sa and we obtained 



Sac 



In {^p\Jk~A + VPc Jk A - 1) 



(3 c Jk A (l - pa) 



In 



1 + Sac 

i-suc 



2Jk A (l -pa)Sac 



(9) 
(10) 



This set of equations determines the critical point 
(T c , Sac) for the first order phase transition between the 
antiparallel and parallel states. 

If we multiply EqJS]by (3 C and Eq[TU]by Sac we get (3 c Sa c 
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FIG. 3: Graphs of Sa(T) for Na = N B = 5000 and k = const 
= 10. The dotted line is for p = (unconnected networks) 
and the rest are for p — 0.1. The dashed line is the graph 
with forced Sa = — Sb (this means forced second order phase 
transition), the solid lines are without such forcing, for par- 
allel and antiparallel initial ordering. The first order phase 
transition is evident for antiparallel case. 
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FIG. 4: Dependence of critical temperature T c on the num- 
ber of inter-network connections Eab for two Barabasi- Albert 
networks. Lines are analytic predictions based on second or- 
der phase transition assumption, while symbols are critical 
temperatures obtained form map iterations. Solid line and 
triangles correspond to Na = Nb = 5000, dashed line and 
squares correspond to Na = 6000, Nb = 4000, while dotted 
line and circles correspond to Na = 8000, TVs = 2000. 



in both and can compare the right sides, obtaining a re- 
lation 

PcJk A + ^p c Jk A (f3 c Jk A - 1) - 1 

ii c = = . (11) 

PcJk A + y/(3 c Jk A ((3 c Jk A - 1) 

Comparing this with Eql9] we obtain a single implicit 
equation for /3 C and p A , that can be simplified to get 

1 + v/1 - \/(J3 c Jk A ) 
Pa = 1 < = ■ 

l3 c Jk A + ^p c Jk A (p c Jk A - 1) - 1 

■]n(y/0 c Jk A + y/0 c Jk A -l) . (12) 

Drawing p A (/3 c ) and changing axes yields a dependence 
of T c on parameter p A (see Figure [5]) . 

We can also approximate the behavior of the solution 
for small p A . Our conditions (Eq|7][E|) can be written 

tanh(/3 c Jk A (l - p A )S Ac ) = S Ac (13) 
cosh 2 (/3 c Jk A (l - PA )S Ac ) = f3 c Jk A (14) 

If we multiply the equations sidewise, we obtain a single 
equation for a multiple X — f3 c Jk A S Ac - 

sinh(2(l-p A )X) =X (15) 

We know that for very small p A the value of S Ac is very 
small, thus X and whole argument of hyperbolic sinus is 
also small and we can approximate it around 

2(1 - p A )X + (2(1 - PA )X) 3 /6 » 2X (16) 

we can calculate the approximate value of X 




Putting the result into the Equation Q3] we obtain follow- 
ing 

cosh 2 ((l- Pyl )V(3/2) PA ) = p c Jk A (18) 

Since the argument of cosh 2 is very small thanks to small 
p A value, we can approximate cosh 2 x = l+x 2 and finally 
obtain T c « k A (l - (3/2)p A ) 

So far, we have concentrated on a case of constant 
node degree k and two networks of same size. Without 
such simplifications, the equations are very hard to solve 
analytically. We have studied more complex cases using 
map iterations and Monte-Carlo simulations. 



V. MAP ITERATIONS 

Since the problem of the first order phase transition 
could not be solved fully analytically, we have used nu- 
merical methods. 

We consider a two-dimensional map 

S A (t + 1) - ]T P{k A )j^ tanh WJA A k A S A (t)+ 

k A 

+f3J BA p A k A S B {t)), (19) 
S B {t + 1) = Y P{k B )-^- tanh (f3J BB k B S B (t)+ 

+PJABP B k B S A (t)) . (20) 

where the S A (t), S B (t) are variables and the rest 
are parameters, including given degree distributions 
P(k A ) and P{k B ). With our definition of weighted 
spin S = l/EJ2ikiSi, where Si are spin values of 
node i, ki are degrees and E is number of edges in 
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FIG. 5: Dependence of critical temperature T c on the param- 
eter p for two constant degree networks k — 10. The straight 
gray line is T c (p) if the transition was of second order. The 
solid line is analytical prediction (Eq|12[l. the plus symbols 
are map iterations, while circles, triangles and diamonds are 
results of numerical Monte-Carlo simulations. Circles are for 
r = 100, triangles up are for r = 30, triangles down are for 
t = 200. Diamonds are for r = 100 but for networks of size 
N = N A = N B = 50000. 



network, it can have values from range [—1,1]. We 
assume Jaa = Jbb = Jab = Jba = J and express 
all temperatures in units of coupling constant J over 
Boltzmann constant k, so we can omit these constants 
in the equations and have j3 — l/T. 
We investigate the dependence of a stable state spin on 
the temperature S(T) in the antiparallel ordering of both 
networks Sa{T) = 1,Sb = — 1- Since the system is fully 
symmetric, below T c we have S(T) = S A (T) = -S B (T). 
At T c , the systems jumps to the parallel ordering. 
Due to deterministic nature of these calculations 
both networks always assume same (negative) spins 
S(T) = S A (T) = S B {T) < 0. By observing S(T) we can 
find the critical temperature T c , where a jump between 
positive and negative spin values occurs (see FigJSJ) . Our 
Sa(T) = (SA{t ma x)) T i S B (T) = {S B {t m ax)) T i where 
tmax = 1000 is the number of iterations of the map 
that have been performed before we assumed it reached 
stationary solution for sure. 

We investigated various T ranges, usually around the 
critical temperature T Cl with the temperature step 
AT = 0.2. Our networks were of size N A = N B = 5000 
and usually possessed a power law degree distribution 
P(k A ) = -P(fcs) taken from a Barabasi-Albert network 
growth simulation or constant degree k — const for 
testing of the analytical equations. Since the networks 
are same pa — Pb = P- 

We have investigated the dependence of the critical 
temperature T c for two networks with constant k = 10. 
The results are in Figure \5\ As can be seen, the 
map iterations do not agree with analytical equations 
exactly. This is probably due to the limited accuracy of 
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FIG. 6: Dependence of |(Sa + Sb)/2| (lines with near zero 
value below T c ) and \Sa\ on the temperature for Monte-Carlo 
simluations. The thick grey lines are for p = 0, solid lines are 
for p = 0.02, dashed for p = 0.1 and dotted for p = 0.2. The 
simulations start from antiparallel ordering Sa = — Sb = 1, 
Na = Nb = 5000, (k A ) = (fcs) = 10 and are performed for 
t = 100. 

numerical calculations, that near such critical point can 
play crucial role. Numerical noise can tip the system 
over the edge into the parallel state. 

We have not however introduced the iterated map to 
investigate what we can analytically. We have performed 
the iterations of the map for the scale-free distribution 
of degrees. The distribution was generated with same al- 
gorithm of Barabasi-Albert network growth as in Monte- 
Carlo simulations, to allow better comparison. 

Looking at Figure [4] it is evident, that the analytic 
results based on assumption of second order phase tran- 
sition are incorrect. For small p, the first order phase 
transition critical temperature is linearly dependent on 
the parameter p, but with different factor. For higher 
inter-network connection number, the dependence is no 
longer linear. 

VI. MONTE-CARLO SIMULATIONS 

Our investigation would not be complete, if we didn't 
use numerical simulations to test our results. We have 
performed Monte-Carlo simulations of the Ising model 
on two inter-connected Barabasi-Albert networks. We 
have simulated networks where N A = Nb = 5000 and 
(k A ) = (kg) = 10. 

The simulation for each temperature T is independent 
on others and starts from an antiparallely ordered 
system S A = —Sb = L The dynamics of Ising model 
are applied for r = 100 time steps and then we perform 
measurement of average \(S A + <Sb)/2| and \Sa\ during 
another r = 100 time steps. One time step equals 
Na + Nb random single node updates, what means 
average one update per node. We have chosen the time 
r so that the network has enough time to relax to the 
equilibrium state, but not enough to reliably jump to 
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FIG. 7: Dependence of weighted spin absolute value \S\ and 
of weighted spin of single network Sa- The results are for 
N A = N B = 5000 and (k A ) = (k B ) = 10. The solid lines 
are for Eab = 0, dashed for Eab = 5000 and dotted for 
Eab = 15000. The upper lines are absolute spins, while the 
bottom, reaching are for single network weighted spin Sa- 
The weighted spin values above 1 result from increasing (k) 
due to interconnections. 



FIG. 8: Dependence of T c on the inter-network link number 
Eab ~ p. The line is analytic prediction of second order phase 
transition model, triangles are map iterations that display 
first order phase transition, while circles are data obtained 
from Monte-Carlo simulations. The plus symbols are same as 
circles, but re-scaled to have same value at Eab — as map 
iterations. 



the parallel state due to the temperature noise. With 
our chosen r = 100 value, the results change little if we 
increase it further, thus we can be sure that the time is 
enough for the network to relax. 

The example of the simulation results are presented 
at Figure El We measure \(Sa + Sb)/%\ and \Sa\ be- 
cause the first order phase transition can be spotted on 
graphs of these values. Since we start from antiparal- 
lel ordering, \(Sa + <Sb)/2| is close to zero below T c , as 
both networks have same sizes and weighter spin S val- 
ues, only of opposite sign, so the total is close to zero. It 
is not exactly zero, only because of fluctuations. Since 
we measure the absolute value, those fluctuations do not 
cancel each other, but add up, resulting in non-zero total 
absolute value of spin. 

When we reach critical temperature T c , 1(5^ + Sb)/2\ 
becomes positive as the networks order parallely. The 
sudden change of total weighted spin means we have 
first order phase transition. We have to use the absolute 
value, since different simulations order either with posi- 
tive or negative spin with equal probability, so if we didn't 
use average of absolute value, we would not be able to see 
the transition point. As the temperature grows higher, 
the networks order parallely with lesser and lesser value 
of spin and finally at temperature T c+ they become para- 
magnetic. This transition is of second order. We do not 
investigate this transition, as it was done before fic| . 
|iSa| behaves similarly, except below T c it is positive. At 
T c there is a sudden change, since the weighted spin of 
networks ordered paralelly is higher at same temperature 
T than the weighted spin of networks ordered antipar- 
alelly. 

We assume that the local minimum of \Sa\ is at T c , just 
before the networks start to order parallely. We inves- 
tigate \(Sa + Sb)/2\ only to confirm that the minimum 



of \Sa\ is indeed at the temperature where the system is 
about to switch to parallel state and \(Sa + 5b)/2| be- 
comes positive. 

We have investigated the dependence of T c on the num- 
ber of inter-network links Eab ~ P- First, we have taken 
the case of k = const, to test how the simulations com- 
pare to analytic results and map iterations. The results 
(Figj5j) indicate, that while critical temperatures T c are 
different than predicted analytically, but the error is not 
large. The fact, that the temperatures drop to zero at 
around p = 0.5, not around p = 1 shows that mean-field 
method does not describe the dynamics of the system ac- 
curately. In real systems, nodes in one network can be 
influenced by the second network stronger than by their 
own. 

Our main results concern the case of the Barabasi- Albert 
networks. The weighted spin against temperature for 
several different interconnection densities is shown in Fig- 
ure [7] 

The jumps of IS"! at values far from zero prove that 
the transition is indeed of the first order as expected, 
since second order phase transition would have weighted 
spin drop to zero (or almost zero, because of fluctua- 
tions) before jumping back to positive (parallel order- 
ing). The dependence of T c on Eab, obtained from data 
that are partially shown at Figure [7] is shown at Figure 
[5] The critical temperature T c is much lower than pre- 
dicted by either analytics or map iterations, but this is 
general problem with Ising model critical temperature in 
B-A network. If for a moment, we omit this and re-scale 
our results to have same value for unconnected networks, 
we obtain relatively good agreement with map iterations. 
Above Eab ~ 15000, the results from map iterations and 
simulations start to differ strongly. This is possibly due 
to simulations having limited system size, so antiparallel 
ordering is quickly destroyed by fluctuations and system 
reverts to parallel ordering that has lower energy. How- 
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ever, for lower interconnection densities, the simulations 
agree with map iterations. 

VII. CONCLUSIONS 

We have shown that in a system of two connected net- 
works, one of two temperature driven phase transitions 
is of first order, unlike classical Ising phase transitions 
that are of second order. The dependence of the criti- 
cal temperature on the interaction strength between the 
networks is complex. The temperatures are lower than a 



theory based on second order phase transition predicts. 
The conclusions are backed up by the numerical simula- 
tions. 
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